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The local dynamical features of a PEO melt studied by MD simulations are compared to two 
model chain systems, namely the well-known Rouse model as well as the semiflexible chain model 
(SFCM) that additionally incorporates chain stiffness. Apart from the analysis of rather general 
quantities such as the mean square displacement (MSD), we present a new statistical method to 
extract the local bead mobility from the simulation data on the basis of the Langevin equation, 
thus providing a complementary approach to the classical Rouse-mode analysis. This allows us to 
check the validity of the Langevin equation and, as a consequence, the Rouse model. Moreover, the 
new method has a broad range of applications for the analysis of the dynamics of more complex 
polymeric systems like comb-branched polymers or polymer blends. 
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I. INTRODUCTION 

Characterizing the complex dynamics of a polymer 
chain in a melt via simplifying models is an important 
problem in polymer physics. For non-entangled chains 
the Rouse model^'- can be regarded as the standard 
model. Within this model a polymer chain is regarded 
to consist of N harmonically bound beads, characteriz- 
ing the intra-molecular forces. Furthermore, the inter- 
molecular interactions are described by stochastic forces, 
uncorrelated in time and space. Considering the over- 
damped limit one can therefore formulate the resulting 
Langevin equation. For bead n of a bead-spring chain it 
is given by 

^^^ = -^ (2R„ (t) - R„+i (t) - R„_i (t)) + iF„ it) , 

(1) 
where C is the monomeric friction coefficient and k — 
3kBT/b^ is the entropic force constant of the springs 
with the average squared length b^. The amplitude 
of the stochastic Brownian force is determined by the 
fluctuation-dissipation theorem^ 



{Fn,a (t) F^^p (0)) = 2kBT(:5n„,8^,t36 [t) 



(2) 



where n and m denote the monomer indices and a and 
P denote the three spatial directions. 

An important observable frequently discussed in the 
context of simulations of polymeric systems and the elu- 
cidation of the quality of the Rouse model is the mean 
square displacement (MSD)1"1. Naturally, the MSD of 
the monomers of a polymer chain shows subdiffusive mo- 
tion on a time scale below the longest relaxation time tr 
due to the chain connectivity. 

Although in general reasonable agreement with the 
very simple Rouse model and simulation data is observed, 



some specific deviations are reported. First, from com- 
parison of the MSD with NSE experiments clear devia- 
tions from Gaussian behavior were reported^'^'^. Second, 
rather than the predicted proportionafity of MSD oc t" 
with a — 0.5 for time scales tr/N"^ < t < tji simula- 
tions of more realistic polymer chains yield a Ri 0.6^"—. 
Third, memory effects for the stochastic forces were 
observed2r22, giving rise to the formulation of more 
compHcated approaches like the generalized Langevin 
equation (CLE) and mode coupling theory (MCT)ii"— . 
Forth, the short-time monomer displacement, which re- 
flects the motion on a local scale, will additionally con- 
tain dynamical features imposed by the detailed chemical 
structure which go beyond the predictions of the Rouse 
model. 

Alternatively, the quality of the Rouse model is probed 
by analysis of the normal modes which in the ideal 
case display a specific dependence on the mode num- 
ber when analyzing their amplitudes and their relaxation 
timegl'— . However, the deviations mentioned above are 
also observed in the behavior of the higher-order Rouse 
modes^i^iiS,. 

The observables, reported so far, are either the normal 
modes themselves or, in case of the MSD, can be ex- 
pressed as a sum over the different normal modes. The 
specific time-dependence of a mode, however, is naturally 
determined by the intra- and inter-chain contributions in 
eq. [U the latter containing the complex interaction with 
adjacent chains. As a consequence, a direct identification 
of the intra-chain contributions is not possible and, as a 
result, the MSD is prone to all different non-idealities. 

More generally, one faces the situation that in order to 
elucidate the quality of the local Langevin equation one 
typically analyzes non-local observables such as normal 
modes or their sum. Here we present a new approach 
(pq-method) which allows us to directly identify the lo- 
cal intra-molecular interactions in eq. [TJ This approach 



is applied to three different model systems. First, we 
analyze an ideal Rouse chain in order to verify our sta- 
tistical approach. Second, we use MD simulations of a 
chemically realistic PEO melt and compare the results 
with the predictions from the Rouse model. Third, as a 
second model system we consider the semiflexible chain 
model (SFCM)iiA?: which is supposed to be superior to 
the Rouse model to describe the dynamics of chemically 
realistic polymers. 

After introduction of the pq-method we mention some 
simulation details for all three systems (PEO, Rouse, 
SFCM). We start with the discussion of standard ob- 
servables, namely the MSD and some characteristics of 
forward-backward correlations. Then we apply the pq- 
method to the three different model systems and charac- 
terize in detail the nature of the intra-chain interaction. 
We conclude by indicating some interesting future appli- 
cations of the pq-method. 



II. PQ-METHOD 

We start by defining the quantities 

R„ {t + At) - R„ (i) 



Pn (i, At) 



At 



and 



q„ (i) = 2R„ (t) ~ R„+i (t) ~ R„_i (t) . 



(3) 



(4) 



Then eq. [T] can be rewritten in the form 

k 1 

^lini^p„ (i, Ai) = --q„ (t) + -F„ 



it, At). (5) 



This is a linear relationship between the vectors p„ and 
q„ with a random component F„. The effective bead 
mobility fc/C can now easily be determined statistically 
by performing a linear regression, i.e. 



A: 



k 

c 



(pq)-(p)(q) (pq) 



(q^) - (q)^ 



(6) 



Strictly speaking one has to consider the limit of small At 
and can average over all inner monomers and all times 
(because the dynamics is stationary). Naturally, for a 
Rouse chain one would recover the effective bead mobility 
as it entered the calculation. 

When applying the Rouse model to a chemically realis- 
tic polymer chain such as PEO one must be more careful. 
As a key problem eq. [T] cannot be applied because of bal- 
listic contributions for very small Ai and non-universal 
local contributions for somewhat longer At. This prob- 
lem can be circumvented as follows. We generalize eq. [5] 
to 



p„ (t. At) = - 



^qn(t) + jF„(t,At) 



(7) 



In this way we introduce a time-dependent friction co- 
efficient C(^^)7 which for small At approaches the bare 



friction coefficient. The slope is determined in analogy 
to eq. mby 



A{At) 



(p(At)q) 



(8) 



Strictly speaking the specific value of p„(i, Ai), re- 
fiecting the dynamics of the n-th monomer during the fi- 
nite time Ai, additionally contains contributions beyond 
its direct neighbors. Exactly these contributions will de- 
crease the relative contribution of q„ (i) to the overall dy- 
namics. This will show up as a decrease of fc/C(Ai) with 
increasing Ai (see below). However, this just refiects an 
intrinsic problem of the Rouse model. The model is ap- 
plicable only for some finite time scale Ai for which the 
Langevin equation in a strict sense is no longer valid. By 
comparison with k/({At) of an ideal Rouse chain one can 
directly check the applicability of the Rouse model to the 
PEO dynamics. 

For a direct comparison of PEO with an ideal Rouse 
chain care must be taken in the definition of eq. |4]because 
Gaussian chain properties have to be assured. Since ad- 
jacent chemical monomers are not independent in their 
orientation, the distance of the two bond vectors entering 
q have to be of the order of the Kuhn length 6k- In order 
to achieve this we rewrite eq. S] 



q„ (i) = 2R„, (i) - R„ 



(i) - R„ 



-Aj^ 



(t) , (9) 



where n' numbers the chemical monomers and An' has 
to be chosen large enough that Gaussian chain properties 
are assured to a good approximation. The value of An' 
can be estimated by the characteristic ratio Coo of the 
simulated chain by making use of the identity 



Co 



&K ^ / ((Rn' — Rn'-A 

"^^ I ((R„. -R„-_i)') 




(10) 



where bo is the average bond length between two chemical 
monomers. For PEO, we obtain An' = 4. The model 
chain described via the SFCM is treated in analogy to 
PEO. Since its bending potential is optimized to refiect 
the structural properties of PEO we again obtain An' = 
4. 



III. SIMULATION DETAILS 

For our analysis, MD simulation data of PEO from a 
previous study have been used^^. A system consisting of 
16 PEO chains with 48 monomers each had been sim- 
ulated in an NVT ensemble with the GROMACS sim- 
ulation package using the two-body effective polarizable 
force field described in refj^. The temperature had been 
maintained at T = 450 K by a Nose-Hoover thermostat. 

An idealized Rouse chain has been simulated via Brow- 
nian Dynamics Simulations. The ideal Rouse chain con- 
sists of A^ = 16 monomers, which corresponds to the 



number of effective beads of the PEO chains with A'' = 48 
monomers and Coo == 3-2 as determined from the MD 
simulations. 

Furthermore, Brownian Dynamics Simulations have 
been also performed for the SFCM. In this model a Rouse 
chain is supplemented by the additional potential 



Ue ({R„}) = ke 



n=2 



Rn) (Rn — R-n-l) 



IR. 



■n+l 



JrVri XVr) 



R„ 



(11) 



Here the value of kg has been chosen such that the char- 
acteristic ratio of the model chain matches that of PEO in 
the MD simulations. Of course, the chains of the SFCM 
contained the same number of monomers {N = 48) as 
the PEO chains. 

Using unit values for the friction coefficient (, the tem- 
perature ksT and the mean squared bond length b^, a 
time-step of At = 0.002 turned out to be sufficiently 
small. 

For the determination of tr, that will be used for nor- 
malization in the following, the autocorrelation function 
of the longest Rouse mode was fitted viai^ 



(X(0)X(i)) = (X2(0))exp(-t/TR,). 



(12) 



In the following calculations, the two outermost beads 
of the Rouse chains and the five outermost monomers of 
the SFCM and the PEO chains were ignored to exclude 
chain-end effects. 



IV. CHARACTERIZATION OF THE 
DYNAMICS 

We start by analyzing the MSD in the c.o.m. -frame 
(fig. [1]). The scaling of the axes guarantees that the 
Rouse curve is independent of the specific choice of the 
parameters &, (', T or N. For later purposes we intro- 
duce Atroi — AdN"^ /tii. Naturally, the subdiffusive dy- 
namics mainly reflects the chain connectivity. Evidently 
the SFCM reflects the PEO dynamics signiflcantly bet- 
ter than the Rouse model, in particular for shorter times. 
However, it is difficult to rationalize the nature of the 
residual deviations of the SFCM for shorter times. A 
similar mismatch of the MSD- related quantity S{q,t), 
i.e. the dynamic structure factor, for the Rouse model, 
the SFCM and an all-atom MD simulation was reported 
in ref ji^. By squaring eq. [7] and expressing the MSD as 
Ai^(p^ (At)) and using (qF) — one obtains the contri- 
bution of the restoring and the effective stochastic forces 
to the overall MSD. The latter is a summation of all 
forces not acting as restoring force within the correlator 
(pq) and therefore also containing the motion of more 
remote monomers and the inter-chain contributions. We 
find for all time scales that the contribution of the effec- 
tive stochastic force is 3 — 10 times higher (depending on 
At) than that of the restoring forces, making the MSD 
very sensitive to the special characteristics of the noise. 
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Figure 1: Mean square displacement (MSD) of the monomers 
in the center-of-mass frame for the Rouse chain, the SFCM 
and the PEO melt. 



A more detailed but related analysis is obtained by 
directly analyzing the properties of the subdiffusive dy- 
namics. For Tpi/iV^ < i < tr, the Rouse theory predicts 
a subdiffusive motion with a bead mean square displace- 
ment proportional to t^'^ . In contrast to this, a propor- 
tionality for the monomer mean square displacement of 
MSD (t) oc i" with a ~ 0.6 is observed in fig. [1] a value 
that is generally found in simulations of more realistic 
polymer models'^"^. This indicates that the backdriving 
forces in a real polymer melt are not only given by the 
chain connectivity, but also infiuenced by chain stiffness 
and intra- and intermolecular excluded volume effects. 
As an investigation tool for this characteristic dynamics 
we define the following quantity 



(ri(^)rU^)) 
(r?(t)) 



(13) 



where ri (t) is the displacement of a given monomer dur- 
ing t and r\ (t) subsequent displacement during the next 
time interval of length t. The MSD (r| (2i)) after 2i can 
also be written as the displacement after two successive 
steps with length (rf (t)) during t 

(r^(2i))=2(r?(0)+2(ri(i)r;(i)), (14) 

where the backdriving force is expressed by the second 
term on the rhs. In cases where the MSD obeys a power 
law. I.e. {rj (t)) = {r'^ (i)) = ct" and (4 (2i)) = c(2i)", 
respectively, the exponent a can be related to the back- 
ward correlation a. By dividing eq. E] by (r^ (i)) and 
rearranging the expression one obtains 

a = l + log2(l + a). (15) 

Figure [5] studies the backward correlation a (left axis) 
and the exponent a (right axis) for the simulated Rouse 
chain, the SFCM and the PEO melt (the dotted line cor- 
responds to a freely diffusing particle). The dynamics of 
the Rouse model for very short times is essentially un- 
correlated to its past, as the beads do not experience the 




Figure 2: Correlation of two successive displacement vectors 
(left axis) and exponent of the MSD calculated by eqn [15] 
(right axis). The dotted line indicates the respective values 
of a freely diffusing particle. 




Figure 3: Scatterplot of the corresponding p- and g- values for 
the PEO chains with An' = 4 and At = 84 ps as well as the 
regression line (eq. [7]) and an anharmonic third-order fit (eq. 

iia. 



connectivity constraint yet (fig. [2|). For larger times one 
observes a long crossover to sublinear diffusion ending in 
a minimum witfi a = 0.5. Note tfiat tlie time window 
for this characteristic motion is rather short due to the 
short chain length {N = 16). In contrast to this, the lo- 
cal bending potential of the SFCM will immediately push 
back a monomer bending too far out of the chain curva- 
ture. Thus, for the SFCM, the motion of a monomer 
is very early dominated by backjumps. For larger time 
scales, a remains nearly constant for the SFCM over sev- 
eral orders of magnitude with slight deviations from ideal 
behavior {a w 0.6 — 0.65). 

To a first approximation the PEO dynamics closely 
resembles the SFCM dynamics in agreement with our 
observations from the MSD. First, the dynamics of the 
PEO monomers crosses over from ballistic to backjump- 
dominated motion at Atroi ~ 0.1 {i.e. approximately 
5ps, during which the monomers move on average 0.75 
times the chemical bond length bo). The minimum at 
Atici ~ 0.1 can be identified as the confinement of 
the PEO monomers in a cage imposed by surrounding 
monomers in the first coordination sphere. The same 
effect is also observed for low-molecular fluidsi^. Sec- 
ond, for Airoi ^ 1 the backward-correlations are again 
somewhat larger. This can at least partly be related to 
the subdiffusivc ccntcr-of-mass motion. For real polymer 
melts, a slightly subdiffusive behavior (with a w 0.8) 
on a time scale shorter than the Rouse time was found 
in simulations^i^iiSii^ and NSE experiment&i^'^ and is 
predicted by theoretical analyses^^'^-^. This is also con- 
firmed in fig. [2] where we also observe a slightly subdiffu- 
sive behavior of the PEO chain's center-of-mass motion 
throughout the entire observation time larger than the 
ballistic regime. The trend of the backward-correlation 
seems to become weaker with increasing time, however, 
we cannot make clear predictions within the error bars. 
The additional backdriving motion acting on the PEO 
chains has been interpreted within the well known cor- 



relation hole picture^^, in which regions of the chain 
that show correlated motion effectively repel each other 
for entropic reasons. From theoretical and numerical 
calculations^ it was shown that the structure of a poly- 
mer chain in the melt can be viewed as a hierarchy of 
correlation holes throughout several length scales, where 
the repulsion is decreasing with increasing length scale. 
Of course, for very long time scales standard diffusive be- 
havior would be observed which is, however, beyond the 
time scales accessible by our simulations. 



V. RESULTS OF THE PQ-METHOD 

As motivated above the pq-method is supposed to yield 
a direct check of the intra-molecular interaction. In par- 
ticular, the specific terms that cause the deviations of the 
MSDs average out in the quantity A (At) . 

Figure [3] shows a scatterplot of the p- and g- values of 
the PEO melt with At — 84 ps as well as the correspond- 
ing regression line with slope A{At) (solid line); see eq. 

To superimpose the curves obtained by eq. [8] we nor- 
mahzed A{At) by the value k/( = N"^ /{■k'^tk) expected 
for the continuous Rouse chain^. As mentioned already 
above. At was normalized by t-h/N'^, characterizing the 
local relaxation time scale within the Rouse model. In 
what follows the normalized time is denoted At^^. 

The normalized curves for ^i.ci(Ati.ci) for a Rouse 
chain, the SFCM the PEO system are shown in fig. |4l 
One observes that from about Atrci ~ 3 (i.e. At « 84 ps) 
the Rouse and the PEO curve agree very well. At this 
time, a PEO monomer has moved approximately twice 
the chemical bond length 6o) showing that at this point 
the microscopic influence of the local potentials is aver- 
aged out and the coarse-grained Rouse picture becomes 
applicable. For shorter At the curve of the Rouse chain 
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Figure 4: Normalized slope A (At) as a function of normalized 
time for the three different systems analyzed in this work. 



converges to the theoretical value of unity. Interestingly, 
around Afrci ~ 3 the value of A{At) for the Rouse sys- 
tem has already reached this theoretical short-time limit. 
One may conclude that at least for PEO the applicabil- 
ity of the Rouse model holds without having to introduce 
renormalized values for the friction coefficient. 

In contrast to the Rouse model, for the SFCM and for 
PEO A (At) further increases with decreasing At. This 
is a consequence of the chain stiffness which enhances the 
effect of a local chain curvature. The PEO curve is in per- 
fect agreement with the SFCM from about Atrei ~ O.lps 
{At « 5ps). Note that this is also the time scale for 
which the forward-backward coefficient a becomes simi- 
lar for both systems. This result clearly shows that the 
dynamical features of a PEO melt are basically captured 
by the SFCM except for shortest times. The deviations 
for At < 5ps are not due to inertia effects because the 
velocity autocorrelation function in the MD simulation 
has already fully decayed around 2ps. Thus, the specific 
properties of the local PEO dynamics is only captured 
by the SFCM for At > 5ps. 

According to the Rouse model, the restoring force act- 
ing on a given bead linearly depends on its elongation (eq. 
[T|) , yielding the characteristic harmonic potential. As can 
already be seen from fig. [Sj an anharmonic fit including 
a third-order term is much more appropriate to describe 
the data (dashed line). In the following, we wish to exam- 
ine the assumption of harmonic restoring forces for the 
SFCM and PEO in more detail. The individual g-values 
were assorted into a histogram, thus yielding the aver- 
ages (q) and (p) of the corresponding p- values for each 
bin (fig. [5]). For sufficiently small bead elongations the 
restoring force indeed is linear. However, when going to 
larger q, this approximation is no longer valid for PEO 
and the SFCM especially in the limit of short At, and 
anharmonicities become noticeable. For the SFCM and 
for PEO, the anharmonicities are similarly pronounced, 
thus again highlighting that the SFCM effectively incor- 
porates many non-idealities of the PEO melt despite its 
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Figure 5: Normalized average p and q values from the his- 
togram of the individual q values for various At. Inset: an- 
harmonicity parameter 7 defined by eq. 1161 



simplicity. The inset in fig. [5] shows the normalized an- 
harmonicity parameter 7 defined by 



p (At) = ~A (At) (q + 7q3 



(16) 



when additionally considering third-order terms. The an- 
harmonicity expressed by 7 for both systems decays on 
the same timescale. At Atj-ci ~ 3, where the average bead 
mobilities of PEO and the Rouse chain become compara- 
ble (fig. |4]), 7 has nearly decayed to zero. Although the 
SFCM and the PEO chains show an anharmonic behav- 
ior at large elongations g, the linear dependence of the 
restoring force as assumed in the Rouse model is valid 
to a good approximation. Deviations become significant 
only on time scales below the Rouse regime. One might 
argue that the nonlinear dependence of p on g leads to 
higher values of A (At) for short At, while for larger At 
and vanishing anharmonicities the curves of PEO and 
the SFCM converge to the Rouse curve (fig. 2]). How- 
ever, we find that the data points for large g-values in 
figs. [3] and [51 where higher-order terms become impor- 
tant, have a statistical weight that is negligible compared 
to the q values near zero. Consistently, the overall value 
for A (At) does not change significantly when the large 
g-values are excluded from the calculation. 



VI. CONCLUSIONS AND OUTLOOK 

In this article, we presented the pq-method to extract 
the effective bead mobility from simulation data of poly- 
meric systems. From a conceptual point of the view the 
main advantage is the direct accessibility of the intra- 
chain interaction effects. In contrast to this, the analy- 
sis of, e.g., the MSD, immediately also involves possibly 
complex inter-chain contributions. In this way it became 
possible to show that the SFCM is a perfect model for 



PEO already for relatively short time scales. Further- 
more we could show that the dynamics of the discrete 
Rouse model is still characterized by the bare friction 
coefficient at the shortest relevant time scale of 84ps. 

However, also from a practical perspective the pq- 
method has a broad range of applications. For example 
it is possible to determine the friction coefficient of the 
n-th monomer. This may be of importance if the friction 
coefficients are heterogeneously distributed. Examples 
are polymer mixtures {e.g. PEO/PMMA) where the lo- 
cal structure of the slow component may determine the 
local mobility of the fast component. More generally, a 
detailed characterization of local dynamic properties be- 
comes accessible. In contrast, observables such as the 
MSD already contain in an uncontrolled way dynamic 
contributions from adjacent monomers or even complete 



Rouse modes, thus invalidating a strictly local analysis. 
Furthermore, the pq-method can be directly applied to 
the analysis of topologically more complex systems {e.g. 
branched polymers) where the local friction coefficient 
may be directly calculated, e.g.^ in dependence of the 
distance to the main chain. 
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